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Preface 


The  original  research  proposal,  from  which  the  work  reported  here  is  an  outcome, 
had  as  its  objectives: 

•  The  general  quantification  and  understanding  of  the  import  for  linear 
elastic  fracture  mechanics  of  three-dimensional  effects  in  crack/surface 
intersection  configurations,  and  in  particular  the  determination  of 
associated  crack  profiles  that  are  critical  in  an  LEFM  sense. 

•  The  development  of  two  independent  numerical  approaches  for  the  attendant 
analysis,  both  of  which  to  be  sufficiently  well  adapted  to  the  class  of 
problems  of  concern  so  as  to  have  high  resolution;  the  two  methods  to  be 
used  jointly  initially  so  as  to  provide  mutual  confirmation  but  ultimately 
being  used  separately  to  each's  respective  greatest  effectiveness. 

The  extent  to  which  the  first  objective  has  been  met  is  sunmarized  in  the 
abstract  that  follows  immediately  and  described  with  greater  detail  and  commen¬ 
tary  in  the  last  of  Section  4  and  Section  5.  The  extent  to  which  the  second 
objective  has  been  fulfilled  is  apparent  from  Section  3,  the  first  part  of  Sec¬ 
tion  4  and  the  last  AFOSR  technical  report  filed.  Both  activities  are  placed  in 
perspective  with  other  related  contributions  to  the  literature  in  the  review  in 
Section  1.  In  sum,  these  various  sources  evidence  a  significant  degree  of  success 
in  achieving  the  basic  aims  of  the  research  program. 
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Abstract 


Several  elastic  configurations  containing  cracks  under  transverse 
tension  which  intersect  a  free  surface  are  investigated.  In  order  to 
ensure  reliable  results  two  independent  numerical  methods  are  employed 
on  a  comparison  problem,  each  method  being  tuned  to  handle  the  special 
features  involved.  The  comparison  provides  confidence  in  other  results 
which  focus  on  the  key  quantity  in  linear  elastic  fracture  mechanics, 
the  energy  release  rate.  These  findings  may  be  summarized  as  follows: 
that  the  decays  in  the  energy  release  rates  found  as  the  free  surface  is 
approached  in  the  various  problems  treated  are  probably  not  significant 
from  a  fracture  toughness  testing  point  of  view  and  not  of  major  consequence 


1.  Introduction 


During  the  past  twenty-five  years  linear  elastic  fracture  mechanics  (LEFM) 
has  emerged  as  the  technology  whereby  the  engineer  predicts  fast  fracture  in 
brittle  instances  and  tracks  cyclic  crack  growth  even  in  some  ductile  materials. 
The  basis  of  this  technology  lies  in  the  identification  of  the  classical 
singular  stresses  and  strains  that  occur  at  a  crack  tip  in  plane  elasticity, 
followed  by  an  integration  of  these  unrealistic  fields  to  form  a  physically 
meaningful,  energy  release  rate,  or  crack  driving  energy,  G:  then  unstable 
crack  propagation  is  predicted  when  G  attains  a  critical  material  value, 
and  fatigue  crack  growth  under  cyclic  loading  can  be  estimated  using  a  data 
reduction  scheme  based  in  essence  on  the  oscillation  in  /G,  a/g.  While  the 
approach  within  its  limitations  is  well  understood  in  two  dimensions  it  is 
less  so  in  three.  Moreover  configurations  occur  quite  often  in  practice 
that  are  inherently  three-dimensional;  for  example,  the  intersection  of  a 
crack  with  a  surface,  especially  in  the  presence  of  shallow  residual  stress 
distributions.  Our  objective  here  is  to  further  the  understanding  of  the 
importance  of  such  three-dimensional  features  within  the  context  of  LEFM. 

Specifically  we  focus  attention  on  cracks  under  transverse  loading  that 
intersect  free  surfaces  because  of  their  potentially  major  implications  in 
engineering,  and  seek  an  appreciation  of  associated  three-dimensional  (3-D) 
effects  on  the  all-important  energy  release  rate.  Thus  we  need  to  identify 

the  local  singular  character  present  in  these  3-D  situations  and,  since  G 
can  no  longer  be  taken  as  being  constant  along  the  crack  front,  to  assess 

Both  of  these  predictive  tools  are  more  usually  couched  in  terms  of  a 
stress  intensity  factor  K  and  its  fluctuation  under  cyclic  loading  AK: 
nonetheless  the  physical  reasoning  underlying  the  approach  remains  that 
of  an  energy  balance  and  K' s  significance  stems  from  its  equivalence  to  G. 


its  variation  there.  In  fact  the  LEFM  postulates  of  a  critical  value  for 
G  for  fast  fracture  and  of  a/g  governing  fatigue  crack  growth  insist  that, 
in  3-D  configurations,  we  look  for  crack-front  profiles  having  a  constant  G 
if  they  are  to  be  everywhere  critical  or  a  constant  a/g  if  they  are  to  be 
associated  with  steady  self-similar  crack  growth.  That  is,  we  need  to  solve 
free  boundary  value  problems  to  determine  the  shapes  of  the  crack  fronts. 

The  consideration  of  these  issues  is  the  major  concern  of  this  paper. 

Over  the  years  a  large  number  of  contributions  to  the  LEFM  literature 
which  address  various  aspects  of  the  Influence  of  three-dimensionality  have 
been  made  -  see  Panasyuk,  Andrejkiv  and  Stadnik  ClD»  C23  for  reviews  through 
1981  which  together  cite  some  500  odd  related  references.  In  commenting 
on  the  impact  of  these  investigations  on  the  issues  of  interest  here,  we 
begin  with  attempts  to  elucidate  local  singular  nature  which  are  primarily 
analytical,  then  turn  to  numerical  treatments  of  more  global  problems. 

To  fix  ideas  in  our  initial  commentary  we  consider  a  generic  configur¬ 
ation  entailing  the  normal  intersection  of  a  free  surface  by  a  crack  under 
transverse  tension  (Fig.  1) .  If  R  denotes  this  cracked  elastic  region  we 
have,  on  introducing  spherical  polar  coordinates  (p,0,$)  with  origin  0  at 
the  crack/surface  intersection  or  crack  vertex, 

R  •  { (p ,©,<}>)  1  0  <  p  <  »,  0  <  6  <  ir/2,  0  <  4>  <  2ir}.  (1) 

To  aid  the  exposition  further  we  also  introduce  cylindrical  polar  coordinates 
(r,+,z)  sharing  the  same  origin  and  related  to  (p , e , 4>)  by 

r  ■  psin8,  ♦  *  ♦,  z  ■  pcos0,  (2) 

on  R.  The  key  question  pertaining  to  this  configuration  is  what  are  the 
possible  singular  aspects  of  the  elastic  stresses  near  the  crack  front  in 
the  vicinity  of  0,  viz.  ,  as  p  +  0  or  as  r,  z  0. 


To  address  this  question  we  first  draw  on  the  work  of  Aksentlan  C3D 
which  expands  the  three-dimensional  equations  of  elasticity  in  local  curvi¬ 
linear  coordinates  to  establish  that  the  stress  behavior  anywhere  along  the 
crack  front  in  R,  except  right  in  the  upper  plane  (8  -  ir/2  or  z  ■  0),  is 
governed  by  the  same  eigenvalue  equation  as  that  derived  by  Williams 
for  the  two-dimensional  case.  Thus  for  this  part  of  the  crack  front  we  have 
the  classical  inverse  square-root  singularity, 

O" -  0( l//x~ )  as  r  +  0  (z  >  0),  (3) 

on  R.  Of  course  there  is  an  infinity  of  other  eigenvalues  giving  rise  to 
eigenfunctions  with  stronger  singularities  than  that  of  (3)  which  are 
acceptable  local  solutions.*  Although  these  singularities  are  usually 
disregarded  on  the  grounds  that  they  lead  to  unbounded  displacements,  the 
dislike  of  such  physically  unrealistic  displacement  responses  does  not 
really  constitute  an  argument  for  their  exclusion  in  any  completely  for¬ 
mulated  global  problem  -  after  all,  the  stresses  admitted  in  (3)  are  cer¬ 
tainly  not  physically  appealing  even  though  their  displacements  are  finite. 

A  real  reason  prohibiting  these  more  singular  stresses  is  that  there  are 
associated  forces  on  plane  regions  of  finite  extent  ahead  of  the  crack  that 
are  infinite ,  so  that,  provided  we  agree  not  to  load  any  finite  surface 
within  R  with  an  infinite  force,  equilibrium  insists  that  they  do  not  par¬ 
ticipate.  Physically,  it  would  be  better  still  if  one  could  argue  against 
the  participation  of  the  classical  singular  stresses  too.  Rather  though, 
the  absence  of  any  two-dimensional  eigenfunction  which  transmits  a  transverse 
tensile  stress  ahead  of  the  crack,  and  does  not  vanish  as  r  -*■  0,  ensures 
that  the  singularity  in  (3)  is  always  excited  in  the  two-dimensional  setting 

^Acceptable  local  solutions  •«  the  .«se  that  they  satisfy  the  field  equations 
and  the  local  homogeneous  bou^  ry  «.  -mditions. 


when  the  crack  is  pulled  apart  and  means  we  can  therefore  expect  its  parti¬ 
cipation  in  any  like  three-dimensional  situation.  Fortunately  it  is  precisely 
this  singular  nature  which  can  result  in  a  positive,  finite,  energy  release 
rate  so  that  the  behavior  in  (3)  is  quite  acceptable  from  an  LEFM  point  of 
view. 

It  remains  to  consider,  for  our  generic  3-D  crack  configuration,  the 
singular  character  actually  in  the  upper  plane.  Three  analyses  which  inves¬ 
tigate  this  are  Benthem  £53»  Kawai,  Fujitani  and  Kumagai  £63,  and  Sinclair 
£73*  All  three  model  their  approach  to  some  degree  on  Williams'  separable 
eigenfunction  analysis  £43,  Benthem  £53  using  spherical  coordinates  and 
satisfying  the  stress-free  conditions  on  the  crack  faces  exactly,  the  stress- 
free  conditions  on  the  upper  surface  approximately;  Kawai  et  at.  £63  using 
the  same  coordinates  but  satisfying  the  crack-face  conditions  approximately, 
the  upper  surface  conditions  exactly;  and  Sinclair  £73  using  cylindrical 
coordinates  and  meeting  the  boundary  conditions  in  a  similar  way  to  Benthem 
£53.  All  three  confirm  the  classical  singular  nature  away  from  the  upper 
surface,  as  in  (3) .  All  three  recover  the  one  known  stress  singularity 
for  a  Poisson's  ratio  v  of  zero  -  the  value  such  that  two-dimensional 
solutions  completely  comply  with  the  local  three-dimensional  requirements  - 
though  Kawai  et  at.  £63  indicate  the  possibility  of  a  stronger  singularity 
as  well  for  this  case.  And,  aside  from  this  last,  all  three  have  qualita¬ 
tively  similar  behavior  as  the  free  surface  is  approached  with 

O'  *  0( p*//r)  as  r  +  0,  (4) 

on  R  in  Benthem  £53»  Kawai  et  at.  £63t  the  exponent  X  varying  from  0  to 
0.17,  0  to  0.32  in  £53t  £63  respectively,  as  v  ranges  from  0  to  0.5,  while 


CT«  0(z  l/r)  as  r  -*•  0,  n  ■  1,  2,..., 


(5) 


on  R  in  Sinclair  C7U»  for  v  ^  0.  No  completeness  argument  for  any  of  the 
above  analyses  is  apparently  available  and  in  fact  their  difference  in  detail 


may  well  be  due  to  a  lack  thereof.*  Nonetheless,  in  sum  it  would  seem 
reasonable  to  expect  behavior  like  that  of  (3),  (4),  (5)  to  be  found  in  any 
3-D  problem  which  becomes  two-dimensional  when  Poisson's  ratio  is  zero,  and 
accordingly  to  expect  a  bounded  energy  release  rate  which  is  positive  except 
right  in  the  upper  free  surface. 

A  related  analytical  treatment  which  contains  the  local  geometry  of 
v  Fig.  1,  yet  attempts  the  ambitious  task  of  the  closed  form  solution  to  a 
truly  three-dimensional,  global,  crack  problem,  is  furnished  by  Folias  in 
C9H  and  has  generated  some  controversy  -  see  discussion  by  Benthem  and  Koiter 
and  author's  closure  ClOH.  The  main  source  of  contention  is  the  appearance 
in  Folias  CO  of  the  following  dominant  singular  behavior  in  the  normal 
stress  component  a  : 


o 

z 


k  cos9  cos(2v+l)9 
p^V  /psinQ 


cos^  +  ...  as  p  +  0, 


(6) 


where  k  is  a  constant  (k  ^  0).  This  singularity  is  stronger  than  (4),  (5) 
and  for  v  >  1/4  gives  rise  to  unbounded  displacements.  As  pointed  out 
earlier,  while  such  displacements  are  perhaps  physically  counter-intuitive, 
they  are  not  in  themselves  a  mathematical  objection  to  the  validity  of  the 
elasticity  solution  in  Folias  C9],  but  rather  raise  the  possibility  that, 
given  the  loading  in  C9D,  the  stress  in  (6)  does  not  meet  equilibrium  require¬ 
ments.  Indeed  inspection  of  (6)  shows  this  to  be  so,  since  equilibrium  in  the 

z-direction  at  the  free  surface  (0  =  ir/2)  requires  3/3z  of  o  to  be  null 

z 

there  (equivalently  3a  /30  to  be  zero),  and  it  is  not  (for  v  t  0) .  Of  course 

z 

*It  might  appear  that  the  two-dimensional  completeness  proof  of  Gregory  C83, 
via  Aksentian's  argument  CO.  applies  here  for  z  >  0;  unfortunately  an  addi¬ 
tional  separation-of-variables  assumption  is  necessary  in  CO  so  that  com¬ 
pleteness  remains  unproven  even  away  from  z  =  0  for  v  i1  0. 


the  method  of  solution  construction  adopted  by  Folias  in  £93  assures  satis¬ 
faction  of  the  equilibrium  requirements  by  his  stress  fields  in  toto,  so 
there  must  be  further  contributions  to  the  stress  field  in  £93,  not  to  date 

explicitly  extracted,  which  combine  with  a  of  (6)  to  restore  this  compliance 

z 

and  in  doing  so  may  even  annihilate  it!  We  conclude  therefore  that  the 
singularity  present  at  the  crack/surface  intersection  in  Folias'  problem  £93 
is  still  to  be  determined  and  consequently  our  expectations  as  to  the  be¬ 
havior  there  remain  unaltered  by  £93- 

The  investigations  £53.  £63.  £73.  £93  do  have  one  thing  in  common  -  they 
all  attest  to  the  difficulty  of  constructing  analytical  solutions  for  3-D 
cracks  and  it  seems  almost  certain  that  the  more  complex  configurations  of 
concern  here  are  intractable  to  purely  analytical  approaches.  As  a  result 
we  look  to  numerical  methods  and  next  review  the  accompanying  literature, 
starting  with  those  contributions  that  are  most  closely  connected  to  the 
aforementioned  analytical  studies. 

The  local  eigenvalue  problem  indicated  in  Fig.'  1,  with  geometry  as  in 
(1)  if  p  <  1  therein,  is  analyzed  in  Benthem  £ll3  via  finite  differences  and 
the  eigenvalues  in  Benthem* s  earlier  paper  £53  confirmed,  the  maximum  dis¬ 
crepancy  between  the  two  determinations  occuring  when  v  =  1/2  and  being  only 
4Z.  Unfortunately  the  numerical  approach  in  £ll3  does  not  readily  admit  to 
more  general  application  and  probably  the  most  obvious  candidate  towards 
this  end  is  the  finite  element  method  (FEM) .  Bazant  and  Estenssoro  £123 
use  a  version  of  the  finite  element  method  on  the  same  local  geometry  as  in 
Benthem  £ll3*  The  most  distinctive  aspect  of  the  FEM  in  Bazant  and 
Estenssoro  £123  is  the  representation  of  the  fields  near  the  crack  vertex 
with  expressions  which  contain  an  adjustable  exponent  governing  the  singu¬ 
larity  there,  the  exponent  being  estimated  by  a  variational 
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principle  -  an  idea  developed  independently  and  somewhat  differently  by  Bazant 
Cl3ll  and  Swedlow  Cl43,  Cl5l.  Using  a  sequence  of  grids,  Bazant  and  Extenssoro 
Cl23  demonstrate  convergence  and,  by  extrapolating  their  results,  also  confirm 
Benthem's  values  for  the  singular  eigenvalue  C53.  As  a  consequence  of  the 
representations  assumed  in  Bazant  and  Estenssoro  Cl23,  this  confirmation  does 
not  constitute  a  completeness  argument  for  Benthem's  forms  in  C53  but  it  is 
quite  reassuring  as  regards  our  expectations  of  singular  behavior,  hence  G, 
in  this  critical  region.  In  addition  in  Bazant  and  Estenssorc  Cl2]  ,  intersec¬ 
tion  other  than  normal  is  considered  and  angles  found  at  which  the  singularity 
exponent  becomes  1/2  (X  =  0  in  (4)),  thus  angles  of  intersection  for  a  straight 
inclined  crack  for  which  positive  energy  release  rates  may  exist  right  in  the 
upper  surface. 

The  success  enjoyed  by  these  local  numerical  analyses  is  not  found  to 
the  same  extent  in  finite  element  treatments  of  more  global  or  complete 
problems  which  contain  crack/surface  intersections,  although  some  worthwhile 
attempts  have  been  recently  reported:  see  Yamamoto  and  Sumi  Cl63»  Yagawa  and 
Nishioka  Cl7],  Cl83  which  use  superposition;  Raju  and  Newman  Cl93,  Sebra 
Pereira  C203,  Tseng  C2l3,  Hilton  and  Kiefer  C223  which  use  special  elements. 

Two  requirements  for  reliable  results  in  this  sort  of  analysis  are  that  the  crack-tip 
representations  of  field  quantities  be  able  to  model  both  singular  and  regular 
behavior  while  retaining  compatibility  with  the  fields  away  from  the  crack- 
tip,  and  that  any  quadrature  scheme  used  should  recognize  the  singular  repre¬ 
sentation  adequately.  To  varying  degrees  the  investigations  in  Cl6]  -  C223 
meet  these  requirements.  However,  all  of  these  papers  have  to  assume  singular 

*Thi8  list  is  not  intended  to  be  comprehensive  but  merely  representative  of 
the  better  finite  element  solutions  for  the  class  of  complete,  elastic,  3D 
crack  problems  considered  here:  we  are  not  aware,  though,  of  any  other 
markedly  superior  FEM  analysis  for  a  problem  of  this  type. 


character  a  priori  and,  presumably  because  of  the  size  of  the  computational 
task  involved  in  global  three-dimensional  geometries,  none  manage  to  carry 
out  the  calculations  on  a  sufficiently  refined  sequence  of  grids  so  as  to 
establish  numerical  convergence  -  an  essential  demonstration  in  view  of  the 
possibility  of  an  inappropriate  assumption  as  to  singular  character,  and  in 
view  of  the  different  predictions  concerning  the  variation  in  the  energy 
release  rate  in  Cl6H  -  C223  with  some  results  showing  an  increase  in  G  in 
the  near  surface,  some  a  decrease,  and  still  others  no  change.  It  would 
seem  therefore  that,  in  the  global  situation,  a  numerical  approach  with  a 
higher  analytical  component  might  be  advisable  if  computational  efficiency 
is  to  attain  levels  that  enable  convergence  to  be  checked. 

A  more  analytical  procedure  is  the  method  of  lines  which  integrates 
displacements  exactly  in  their  own  directions  and  uses  finite  difference 
approximations  for  their  gradients  in  transverse  directions.  Application 
of  the  method  to  global  configurations  containing  cracks  with  normal  inter¬ 
sections  oi  free  surfaces  is  described  in  Gyekenysei  and  Mendelson  [23]. 

In  C23U  G  can  be  inferred  to  remain  constant  or  to  increase  as  the  free 
surface  is  approached  depending  on  overall  specimen  geometry,  clearly  a 
contradiction  to  the  expectation  expressed  in  (4) ,  (5) .  What  is  not  clear 
in  Gyekenysei  and  Mendelson  C  23U  is  how  well  the  singularities  in  the  dis¬ 
placement  derivatives  are  approximated  by  the  finite  differences,  so  that 
there  is  some  doubt  as  to  the  reliability  of  the  results  obtained  by  means 
of  this  method  at  this  time. 

Another  more  analytical  approach  than  FEM  is  the  boundary  integral 
equation  method  (BIE)  which  employs  a  Green's  function  to  meet  the  field 
equations  exactly,  thus  reducing  the  problem  to  boundary  condition  satis¬ 
faction  on  surfaces;  i.e.,  the  region  to  be  discretized  for  numerics  is 
reduced  from  three  dimensions  in  FEM  to  two  dimensions  in  BIE.  An  early 


use  of  this  method  on  a  complete  3-D  crack  problem  is  that  by  Cruse  and 
Van  Buren  £24]  and  remains  even  today  one  of  the  more  refined  analyses  -  at 
least  as  judged  by  the  size  of  the  discretization  regions  near  the  crack. 

A  more  recent  application  of  BIE  to  a  problem  of  this  class  is  that  of  Tan 
and  Fenner  £25]  which  has  a  higher  order  representation  of  the  unknowns  in¬ 
volved  yet  still  does  not  take  their  singular  behavior  into  account.  Both 
Cruse  and  Van  Buren  £24]  and  Tan  and  Fenner  £25]  find  evidence  of  a  decay  in 
the  energy  release  rate  as  the  free  surface  is  approached  but  neither  fully 
confirms  their  findings  with  a  convergence  check  on  a  discretization  sequence. 
In  order  to  do  this  check,  special  integral  equations  that  exploit  a  match 
between  a  particular  Green’s  function  and  a  restricted  class  of  three- 
dimensional  crack  geometries  may  prove  advantageous.  Recently  two  such  in¬ 
tegral  equation  sets  have  been  derived  by  Folias  [26]  and  Smetanin  and  Sobol 
£27];  currently  though  no  results  from  their  numerical  analysis  seem  to  be 
available. 

In  the  light  of  the  preceding,  this  investigation  uses  a  special  set 
of  integral  equations  developed  in  Sinclair  £28]  for  an  elastic  half-space 
weakened  by  a  crack  of  finite  width  and  infinite  depth  subjected  to  trans¬ 
verse  excitation.  This  integral  Equation  set  has  a  single  physical  unknown, 
is  more  compact  than  that  of  Folias  £26]  but  comparable  in  complexity  to 
that  of  Smetanin  and  Sobol  £27],  and  has  the  attribute  of  recovering  exactly 
the  closed  form  two-dimensional  solution  that  obtains  for  the  configuration 
under  the  simplifying  assumptions  of  normal  intersection,  uniform  far-field 
tension  and  a  Poisson's  ratio  of  zero.  Furthermore,  the  singular  nature  of  the 
unknown  in  the  set  needs  only  to  be  dominated,  not  estimated  precisely,  to 
ensure  convergence,  and  this  last  concern  is  examined  on  a  sequence  of  three 
discretizations,  the  finest  of  which  provides  information  at  distances  in 
advance  of  the  crack  front  which  are  about  a  factor  of  four  closer  than 
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the  closest  in  Cruse  and  Van  Buren  C24],  On  the  other  hand,  little  is 
learned  of  the  precise  singular  behavior  in  the  crack  vertex  proximity  from 
the  system.  To  compensate  we  also  employ  a  finite  element  method  to  analyse 
a  finite  geometry  of  sufficiently  gross  dimensions  so  as  to  simulate  the 
half-space  implicitly  required  in  our  integral  equation  development.  The 
FEM  features  an  adjustable  singularity  exponent  at  the  crack  vertex  (after 
Bazant  and  Estensscro  Cl2],  Swedlow  Cl5l)  so  that  the  singularity  there  is 
not  completely  assumed  in  advance;  further  the  special  crack-tip  elements 
involved  also  contain  regular  fields  and  are  compatible  with  host  elements, 
and  singular  quadrature  is  used.  By  comparison  with  the  integral  equation 
results  its  convergence  can  be  Inferred.  With  the  two  numerical  approaches 
thus  verified  they  can  then  be  used  to  each's  best  advantage  to  study 
critical  crack  profiles,  residual  stress  effects,  and  the  influence  of 
plate  thickness. 

We  begin  in  Section  2  by  stating  the  3-D  problems  for  analysis  via  the 
Integral  equation  and  finite  element  approaches.  In  Section  3  we  describe 
the  analysis,  commencing  with  the  derivation  and  numerics  associated  with 
the  integral  equations,  and  closing  with  the  development  of  the  FEM.  In 
Section  4  we  present  results  obtained,  commencing  with  those  which  serve 
to  validate  the  overall  treatment,  and  closing  with  those  which  quantify 
selected  three-dimensional  effects  in  LEFM.  Finally  in  Section  5  we  offer 
some  concluding  remarks. 

2.  Formulation  of  three-dimensional  problems 

Here  we  consider  two  companion  geometries  for  the  examination  of  the  three- 
dimensional  aspects  of  crack/surface  interactions.  The  first  is  an  infinite 
geometry  chosen  with  a  view  to  facilitating  integral  equation  development, 
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the  second  a  finite  geometry  designed  to  enable  comparison  of  an  FEM  treat¬ 
ment.  We  start  with  a  formal  statement  of  the  simplest  problem  analyzed  via 
integral  equations  and  continue  with  a  discussion  of  some  generalizations  of 
this  configuration  to  address  other  issues  of  physical  importance.  We  then 
follow  the  same  pattern  in  describing  the  problems  for  finite  element  analysis 

One  configuration  amenable  to  integral  equation  development  is  that  of  an 
elastic  half-space,  containing  a  crack  of  finite  width  and  infinite  depth, 
under  uniform  transverse  tension  at  infinity  o°  (Fig.  2).  This  choice  is 
not  uniquely  well  suited  to  our  present  purpose  but  is  sensible  in  that  it 
has  a  characteristic  length,  is  subjected  to  opening  (Mode  I)  excitation, 
and  does  enable  a  compact  set  of  integral  equations  to  be  subsequently  de¬ 
rived  in  the  next  section.  With  the  development  of  the  integral  equation 
system  in  mind  we  subtract  a  uniform  transverse  tension  of  magnitude  o® 
throughout  an  uncracked  half-space  to  recover  a  configuration  wherein  the 
crack  is  opened  by  a  constant  normal  pressure  equal  to  o°,  then  exploit  the 
symmetry  about  the  crack  plane  to  formulate  the  problem  in  rectangular 
cartesian  coordinates  (x,y,z)  on  the  quarter  space 

Q.  ■  {(x,y,z)|  -  «•  <  x  <  »,  0<y<»,  0  <  z  <  «}.  (7) 


For  this  quarter-space  we  denote  the  upper  surface  by  3^0.  and  the  surfaces 
corresponding  to  the  crack  face  and  the  connected  material  in  the  crack 
plane  by  Q.t  3^0.  respectively.  Thus 

3 jQ.  *  { (x,y,z)  |  -  »  <  x  <  <*>,  0  <  y  <  ®,  z  =  0}, 

32£  "  { (x,y, z) |  0  <  |x|  <  a,  y  =  0,  0  <  z  <  »},  (8) 

33 £  *  { (x,y, z) |  a  <  |x|  <  ®,  y  =  0,  0  <  z  <  »}, 


wherein  a  is  the  semi-crack  length.  We  seek  then,  in  general,  the  stresses 


0  y  0  *  0  * 

x  y  z 


x  ,  x  and  displacements  u,  v,  w  as  functions  of  x,  y,  z 
yz  zr 


x 


throughout  Q.,  satisfying:  the  three-dimensional  stress  equations  of 
equilibrium  in  the  absence  of  body  forces,  on  Q,;  the  stress-displacement 
relations  for  a  homogeneous,  isotropic  and  linear  elastic  continuum  with 
shear  modulus  y  and  Poisson's  ratio  v,  throughout  £;*  the  stress-free  condi¬ 
tions  on  the  upper  surface, 

a  m  r  *  t  =0,  (9) 

zyzzx 

on  djQ;  the  opening-pressure  stress  conditions  on  the  crack  face, 

a  "  — o°  ,  t  =  t  =0,  (10) 

y  y  xy  yz 

on  the  symmetry  conditions  ahead  of  the  crack, 

t  =  t  =  0,  v  «  0,  (11) 

xy  yz 

on  and  finally  the  conditions  at  infinity  which  insist  that  the  stresses 

approach  a  plane  strain  state  for  a  crack  opened  under  constant  pressure, 

_  _  / 1  \  _  _  ./1\  _  _  ./1\ 2  .  __ 


O(l),  T 

■  0(1),  T  »  i 

xy 

xz 

0(1.),  T 

=  0(1),  T= 

xy 

yz 

aP  +  0(1), 

T  =  0(1),  T 

yz 

on  where  ap  is  the  plane  strain,  out-of-plane,  normal  stress  (refer,  for 
z 

instance.  Green  and  Zerna  C303,  §8.17  for  details), 

x  $  <j> 

ap  -  2v  o°  l  r=  cos  13,  (13) 

z  y  /rjr2  °  2  2 


r±  -  *t(x  +  i(5-3i)a/2)2+y2],  =  sin_1(y/ri) (i  =  0,1,2), 

on  In  particular  we  wish  to  draw  from  the  stresses  and  displacements 
complying  with  this  formulation  the  attendant  energy  release  rate  and  con¬ 
sider  three-dimensional  effects  thereon. 


*The  usual  notation  for  stresses  and  displacements  is  employed;  for  details 
of  the  field  equations  see,  for  example,  Timoshenko  and  Goodier  C293,  pp.  7, 
,236. 
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The  energy  release  rate  G  here  may  be  represented  using  any  one  of 
several  equivalent  integral  forms.  Viewing  G  as  the  incipient  local  work 
required  to  heal  the  direct  extension  of  a  crack  in  the  opening  mode  per 
unit  length  of  that  extension,  one  such  integral  is,  after  Irwin  C3lD, 


G 


11m 
6a  -*■  0 


a+6a 

—  I  oa  6vdx. 
6a  J  y 


(14) 


Here  6 a  is  the  extension  of  the  crack  in  the  x-direction  through  material 
initially  experiencing  the  tensile  stress  a and  6v  is  the  opening  dis¬ 
placement  accompanying  this  extension.  Implicit  in  (14)  and  its  equivalent 
forms  is  a  requirement  on  the  local  singular  character  so  that  the  integra¬ 
tion  yields  a  finite  non-trivial  result.  The  determination  of  this  requisite 
singular  nature  proceeds  routinely.  We  take  the  most  singular  aspect  of  the 

d 

stress  distribution  o  along  the  x-axis  to  be  described  by, 

y  * 

oa  *  C1  (x-a)*  (1  +  0(1))  as  x  -*■  a+  (X  >  -1),  (15) 

where  X  is  the  exponent  characterizing  the  singularity  and  a  constant. 

This  in  turn  implies  an  associated  shifted  displacement  of 

6v  **  C2(a+6a-x)*+^  (1  +  0(1))  as  x  (a  +  6a)  ,  (16) 

C2  a  constant.  Then  introducing  the  forms  (15) ,  (16)  into  (14)  and  making 
the  change  to  the  variable  x'  »  (x-a)/<$a  gives 

G  "  SaTo  C1  C2  6a2X+1B(*+l,*+2),  (17) 

where  B(X+l,X+2)  is  the  Beta  function,  bounded  for  X  >  -1.  Equation  (17) 
insists  that  X  *  -h  in  order  for  G  to  be  bounded  and  positive. 

Also  implicit  in  (14)  is  that  in  two  dimensions  G  is  uniform  along  the 
crack  front,  whereas  in  three  dimensions  this  integral  represents  an  energy 
release  rate  per  unit  length  of  crack  front  and  its  constancy  is  not 
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assured.  Thus  in  order  that  the  entire  crack  front  be  critical,  we  cannot 

assume  its  shape  a  priori  but  must  determine  it  such  that  G  is  uniform. 

Currently  two-dimensional  testing  assumes,  in  effect,  a  straight  crack 

front  to  be  critical  and  on  this  basis  estimates  the  critical  value  G  for 

c 

fracture,  the  material  property  one  seeks  essentially  in  fracture  toughness 
testing.  On  the  one  hand  the  attendant  curvature  restrictions  in  these  tests 
may  result  in  unnecessarily  discarding  useful  data  (as  in  Kaufman  [323) ;  on 
the  other,  the  average  G  found  via  two-dimensional  analysis  for  a  straight 
or  nearly  straight  crack,  set  equal  to  G£  at  fracture,  may  be  significantly 
lower  than  the  actual  three-dimensional  maximum,  thereby  leading  to  an  overly 
conservative  estimate  of  this  critical  energy.  Now  the  usual  physical  sit¬ 
uation  preceding  fracture  -  cyclic  growth  up  to  the  point  of  unstable  crack 
propagation  -  indicates  near  uniform  advances  of  curved  crack  profiles  (as 
in,  for  example,  Bell  and  Feeney  [333).  Thus  a  more  realistic  calculation 
of  Gc  would  result  if  we  considered  these  fronts,  which  may  reasonably  be 
taken  to  possess  constant  energy  release  rates,  and  at  failure,  a  constant 
Gc*  In  order  to  state  the  associated  free-boundary  problem  precisely,  we 
need  to  choose  the  plane  in  which  the  energy  release  rate  is  calculated  to 
best  model  the  observed,  approximately,  self-similar,  crack  growth.  Here 
two  likely  candidates  present  themselves;  normal  to  the  crack  front  and 
parallel  to  the  free  surface  (this  last  actually  being  the  direction  in  (14)). 
To  select  one  of  these  we  consider  the  consequences  of  local  crack  advances 
governed  by  a  Paris  data  reduction  scheme  for  cyclic  loading.  For  either 
direction,  beginning  with  an  assumed  profile  of  constant  G,  the  crack  under 
a  Paris  relation  is  advanced  a  uniform  increment  6a  in  a  single  cycle  directed 
in  the  plane  in  which  G  is  calculated.  We  see  then  (Fig.  3)  that  only  a 
crack  advance  associated  with  G  calculated  in  the  plane  parallel  to  the  free 
surface  produces  a  crack  that  retains  its  shape.  Accordingly,  we  choose 
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Figure  3.  Candidate  directions  of  crack  advance  to  achieve  self-similar  profiles. 


this  plane  for  the  calculation  of  the  energy  release  rate.  As  a  result,  the 
pertinent  free  boundary  problem  seeks,  in  addition  to  the  stresses  and  dis¬ 
placements  satisfying  all  of  our  preceding  requirements,  to  determine  the 
crack  profile,  a  =  a(z),  such  that,  in  light  of  the  singularity  requirement 
demanded  by  (17) ,  we  have, 

lim  a  /(  x  -  a)  =  C  ,  (18) 

4>  V  o 

x  -►  a  J 


on  3,Q.,  C  being  a  constant. 
J  o 


Another  distinctly  three-dimensional  problem  in  LEFM  that  can  be  addressed 
through  a  generalization  of  the  previous  integral  equation  formulation  is  the 
effect  on  the  energy  release  rate  of  a  crack  embedded  in  a  non-uniform  residual 
stress  field.  In  practice,  one  attempts  to  induce  near-surface,  compressive, 
residual  stress  distributions  to  inhibit  fatigue  crack  initiation  and  propa¬ 
gation.  It  is  not  intended  here  to  trace  the  actual  physical  processes  which 
lead  to  these  distributions,  a  difficult  non-linear  elasto-plastic  problem, 
but  instead  we  look  to  perform  a  series  of  numerical  experiments  to  quantify 
the  LEFM  implications  of  some  presumed,  representative,  compressive  residual 
stresses.  Based  on  typical  experimental  results  for  compressive  stresses 
imposed  by  a  shot  peening  (see,  for  example,  Almen  and  Black  C34],  p.  58), 
we  choose  to  represent  the  residual  stress  field  as  exponentially  decaying 
away  from  the  free  surface.  More  precisely,  we  take  the  residual  stress 
distribution  on  the  crack  face  to  be 


-C  e 
r 


-yz/t 


(19) 


on  32Q.»  where  the  magnitude  and  depth  of  penetration  of  this  field  are  governed 


by  the  positive  constants  C and  y  respectively.  The  associated  problem  is 


as  previously  stated  except  that  the  crack  is  now  opened  by  the  constant 


pressure  o®  minus  a i.e.,  (10)  becomes 

o  "  -0°  -  aT ,  x  =  t  =  0,  (20) 

y  y  y’  xy  yz 

on  320.*  In  an  attempt  to  ensure  no  interpenetration  of  the  opposing  crack 
flanks  we  make  nowhere  positive  on  by  setting  Cr  equal  to  a®  ,  and 
examine  the  effects  of  different  decay  rates  y. 

Turning  to  problems  to  be  analyzed  via  the  FEM  approach,  the  primary 
one  for  comparison  requires  that  we  select  a  finite  geometry  to  simulate 
the  baseline  integral  equation  configuration  of  Fig.  2 .  Thus  we  make  the 
finite  problem  large  enough  so  that  the  effects  of  its  limited  extent  are 
negligible  at  the  crack  front.  We  arrive  at  the  in-plane  dimensions  by  con 
sidering  the  two-dimensional  analysis  of  Isida  Q35]  for  a  center-cracked 
plate  in  the  opening  mode.  In  C35l  a  crack  length  to  plate  width  ratio, 
a/b,  of  1/10  and  a  plate  length  to  width  ratio,  c/b,  of  2  yields  a  dif¬ 
ference  in  the  singular  participation  factor  of  less  than  one  percent  from 
the  corresponding  two-dimensional  infinite  plate  problem.  Hence  we  take 
these  to  be  the  dimensions  here.  For  the  out-of-plane  dimension,  no  two- 
dimensional  problem  readily  enables  a  similar  sizing.  We  choose  a  plate 
thickness  to  the  crack  length  ratio,  h/a,  of  2  in  an  attempt  to  recover  the 
generalized  plane  strain  solution  at  the  plate  midplane;  our  selection 
requires  subsequent  checking  to  see  if  this  is  in  fact  the  case.  The  re¬ 
sulting  geometry  is  subjected  to  a  uniform  uniaxial  tensile  load  o®  normal 
to  the  crack  plane  on  the  plate  ends  (Fig.  4,  not  to  scale). 

Looking  to  minimize  the  region  to  be  discretized  ultimately  in  the 


finite  element  analysis,  we  exploit  the  symmetries  about  the  plate  mid-face 
the  crack  plane  and  the  plane  normal  to  and  bisecting  the  crack  plane, 
to  formulate  the  baseline  problem  in  rectangular  coordinates  (x,y,z)  on 
the  plate  octant  P, 


Figure  4.  Coordinates  for  the  basic  finite  configuration 


P  ■  {(x,y,z)|  0<x<b,  0  <  y  <  c,  0<z<  h}. 


(21) 


having  plate  upper  face  3^P  and  mid-face  82?  given  by 

3jP  =  {(x,y,z)|  0<x<b,  0<y<c,  z  =  (i-l)h)  (i  =  1,2),  (22) 

with  mid-plate  and  outer  plate  edges  S^P,  9^P  respectively  as  in 

a^P  =  {(x,y,z)j  x  =  (i-3) b,  0<y<c,  0<z<h}(i=  3,4),  (23) 


and  having  crack  flank  S^P,  crack  plane  connecting  surface  3^P,  and  plate 
end  3^P  defined  by 

3±P  =  (  (x,y,  z)  |  (i-5)  (7-i)  a  <  x  <  (6-i)  (7-i)-^~^-  +  b. 


(24) 

y  «  (i-6) (i-5)|,  0  <  z  <  h)  (i  =  5,6,7). 


We  seek  then,  in  general,  the  stresses  and  displacements  as  functions  of 
x,y,z  throughout  P  satisfying:  the  same  elastic  field  equations  as  previously 
outlined,  throughout  P;  the  stress  boundary  conditions  on  the  loaded 
end  and  the  free  crack  flank. 


o  «=o°,x  =  t  *0  and  0  -  t  =  t  =0, 

y  y  xy  yz  y  xy  yz 

on  3?P,  35P  respectively;  the  free  edge  and  free  face  conditions, 


(25) 


O  =  T 

x  xy 


x  =0  and  0  =  t  =  t  =0, 

zx  z  zx  yz 


(26) 


on  34P,  3^P  respectively;  the  symmetry  conditions  on  the  plate  mid- face 
and  the  mid-plate  edge , 


Tyz  -  Tzx  -  0,  V  -  0  and  Txy  0, 


(27) 


on  32P*  a3P  respectively;  and  the  symmetry  conditions  on  the  plane  ahead  of 
the  crack , 


xy 


yz 


0,  v  -  0, 


(28) 


on  a^P.  In  particular  we  wish  to  extract  from  the  stresses  and  displacements 
complying  with  the  above,  the  energy  release  rate  G  as  defined  in  (14). 


In  the  selection  of  the  finite  geometry  we  chose  the  thickness  ratio 
h/a  so  as  to  allow  a  comparison  with  the  basic  infinite  problem  posed.  We 
cannot  make  this  choice  in  advance  but  must  consider  a  succession  of 
problems,  where  we  vary  the  h/a,  until  we  obtain  the  desired  midplane  result. 
Then,  in  addition  to  determining  a  useful  geometry  for  comparison,  we  may 
examine  the  effects  of  plate  thickness  on  the  crack  driving  energy  along  the 
crack  front,  that  is,  the  interaction  of  two  crack/surface  intersections.  With 
the  variation  in  G  quantified,  a  natural  extension  of  the  FEM  investigation 
is  a  search  for  critical  profiles  with  a  constant  G,  to  wit,  profiles  meeting 
(18).  This  last  completes  the  set  of  problems  we  consider. 

3.  Integral  equation  and  finite  element  analyses 

In  this  section  we  establish  a  set  of  integral  equations  for  our  basic  infinite 

3-D  problem  and  discuss  extensions  of  the  set  to  accommodate  the  generaliza¬ 

tions  of  interest.  Then  we  examine  the  issues  entailed  in  the  numerical 
solution  of  such  integral  equation  sets.  Finally  we  describe  the  finite 
element  analysis  of  the  companion  finite  3-D  problems,  paying  particular 
attention  to  the  special  features  involved  in  their  treatment. 

To  improve  the  convergence  of  the  integrals  in  our  integral  equations, 

we  start  by  modifying  the  basic  infinite  problem  by  removing  the  plane  strain 

A 

solution  for  a  crack  opened  by  a  constant  pressure  o°.  Then  the  solution 
to  the  outstanding  three-dimensional  problem  in  essence  becomes  a  local 
free-surface  correction  to  the  two-dimensional  problem.  Formally  the  basic 
problem  statement  for  the  infinite  geometry  remains  intact  except  that:  the 
upper  surface  stress-free  conditions  (9)  and  the  crack-face  opening-pressure 

* 

Indeed  for  the  integral  equation  set  ultimately  constructed  here,  this 
removal  is  essential  if  the  improper  integrals  involved  are  to  converge. 
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conditions  (10)  become 

o  -  -  oP,  t  ■  x  =  0,  (29) 

z  z  yz  xz 

on  (oP  as  in  (13))  and 
l  z 

o  -  t  *=  x  *  0,  (30) 

y  xy  yz  v  ' 

on  82^  respectively,  and  the  infinity  conditions  (12)  are  simplified  with 

all  the  stresses  now  being  0(1). 

To  derive  integral  equations  for  this  modified  problem  we  seek  to 
superimpose  point  loads  on  the  quarter-space  Q [.  Unfortunately  no  solutions 
expressed  in  elementary  closed  forms  appear  to  exist  for  such  point  loads. 
Fortunately  though,  if  we  form  the  quarter-space  as  the  intersection  of  two 
half-spaces  (Fig.  5),  we  have  at  our  disposal  the  Boussinesq  solution  for  a 
normal  point  load  on  a  half-space  (the  simple  forms  for  which  can  be  found  in 
Timoshenko  and  Goodier  £29],  pp.  398-402).  That  is,  we  form  as 

a  -  H,n  h2, 

H.  -  {(x,y,z)|  —  <  x  <  <*»,  -»  <(i-l)y  +  (2-i)z  <  °°, 

(31) 

0  <  (2-i)y  +  (i-l)z  <  «>}  (1-1,2), 
with  corresponding  surfaces  8-  H  given  by 

3 .H  «  {(x,y,z)|  -00  <  x  <  °°,  -°°  <  (i-l)y  +  (2-i)z  <  00 , 

(32) 

(2-i)y  +  (i-l)z  -  0}  (i-l,2). 


Then  the  symmetric  superposition  of  four  equal  point  loads  P  at  x  *  ±a, 

I 

z  -  ±8  on  together  with  four  equal  point  loads  P  at  x  -  ±x,  y  -  ±\l> 

on  <*2^  satisfies  the  field  equations  throughout  £  and  the  zero  shear  stress 

conditions  on  3jQ,  820.,  snd  reflects  the  last  plane  of  symmetry,  the 

yz-plane:  it  remains  to  distribute  P,P  so  as  to  satisfy  a  =  -aK  on  3.Q., 

z  z  1 

a  -  0  on  3 _(£,  and  v-  0  on  3_Q,.  We  now  treat  each  of  these  in  turn. 
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By  superposition,  the  normal  stress  a induced  on  the  symmetry  plane  of 
Hj,  HjOa  2 W,  in  response  to  four  Boussinesq  point  loads  of  magnitude  P 
symmetrically  exerted  at  x  ■  ±a,  2  »  ±8  on  9  is 


o  -  PF(a» $;x,y) , 
z 


on  3  {£  »  W1D32^»  where 


A  l  r  ,  5 


F(a,8;x,y)  =  -  Z  C-  38  y  /  + 

i=l  1 


(l-2v)((y2  +  82)(y  +  2R1)  -  R*)/(rJ  (y  +  R±)2)](,34) 

i  \2  ,  „2  ,  2-i  >i  _  1 


\  =  /[(x  +  (-)  a)z  +  8  +  y  j  (i  =  1,2). 


Hence  if  P  is  replaced  by  the  loadette  p(a,8)  acting  on  the  element  dad6 
centered  on  (a,8X  and  p(a,8)  distributed  throughout  the  quadrant 
(0  <  a  <  »,  0  <  8  <  ■»)  we  have 


p(a*8)F(a, 8;x,y)dad$, 


on  Since  becomes  H ^  on  interchanging  y  and  z,  the  stress  oy  induced 

on  the  symmetry  plane  at  O  in  response  to  four  point  loads  of 

t 

magnitude  P  at  x  =*  ±x,y  =  on  82^,15  given  by  (33)  with  y  interchanged 

f 

with  z,  P  replaced  by  P  ,  and  a,g  by  x.'l'-  It  follows  then  that  the  response 

to  loading  p  (x.'P)  distributed  throughout  (0  <  x  <  00  ,  0  <  i|>  <  ®)  is: 

00 

’  ' 

o  -  p  *(X  ,t|»;x,z)dxdij;,  (36) 


°n  “  «2d  hH-  Now  combining  t^e  distribution  p  on  3 with  the 

response  to  p  there  (35),  satisfaction  of  the  first  of  the  outstanding 
requirements  implies  that 


°z  “  ~P  (x,y)  +  p(a,8)F(a,8;x,y)dad6  **  -  aP, 

z 


0 


on  323. 

Turning  to  the  last  of  the  outstanding  conditions  to  be  met,  the  dis¬ 
placement  condition,  we  observe  that  only  the  distributed  loading  p  on  3^H 
contributes  to  the  prescribed  displacement  component.  To  get  this  contribu¬ 
tion  we  draw  on  the  Boussinesq  solution  to  assemble  the  displacement  in  the  y-dir- 
ection  on  the  plane  1°  response  to  a  set  of  four  equal  point  loads  at 

x  ■  ±a,  z  ■  ±3.  Then  the  requirement  that  the  response  to  the  distribution 
p  be  zero  yields  the  displacement  integral  equation ,  our  second  integral 
equation  in  the  unknown  p. 


p(a,B)f  (a,B;x,z)dctdg  «  0, 


on  3 where 


2  2 

f(a,$;x,z)  =  ~~  l  l  1//  C(x  +  (-^a)2  +  (z  +  (-)J3)2]. 
2mJ  i-1  j=l 


(40) 


(41) 


Equations  (39),  (40)  constitute  our  fundamental  set  of  Fredholm  integral 
equations  in  the  single  unknown  pressure  p;  in  essence,  they  provide  the 
local  three-dimensional  free-surface  correction  to  the  two-dimensional  problem. 
Before  moving  on  to  their  numerical  solution,  some  comments  on  certain 
analytical  aspects  are  in  order. 


Observe  that  oK  is  0(1/ (x  +  y  ))  on  3,0  as  (x  +  y  )  -*■  <*>  and,  since  the 
z  x 

2  2  2  2  2  2 

decay  of  the  Boussinesq  stress-kernel  is  0(1/ (x  +  y  +  z  ))  as  (x  +y  +  z  ) 

2  2 

-*•  the  decay  of  the  unknown  p  can  be  inferred  as  0(1/ (x  +  y  ))  on  32^.  U  3^0. 
2  2 

as  (x  +  y  )  -*■  «•;  under  these  conditions  the  improper  integrals  in  equations 
(39), (40)  are  guaranteed  to  be  convergent.  Note  also  that,  when  v  =  0,  of 

(13)  is  zero,  which  leads  to  p  being  zero  in  (39),  (40)  and  the  only  known 
closed-form  3-D  solution  to  (9) -(12)  is  recovered  exactly  -  the  instance  where 
the  three-dimensional  solution  becomes  two-dimensional.  Unfortunately  equa¬ 
tions  (30),  (40)  appear  resistant  to  the  generation  of  further  complete 
analytical  solutions.  Nonetheless,  by  making  successive  approximations  on 
the  unknown  p  *  p^n\  commencing  with  the  trivial  p^°^  =  0,  and  changing 
variables  by  setting  x-a  *  psin0  and  z  *  pcos0,  p  being  the  distance  from 
the  crack  vertex  and  0  the  angle  away  from  the  crack  front,  we  obtain 


Here  0  is  a  function  of  0  alone  ,  not  everywhere  zero.  Of  course  this 
asymptotic  treatment  is  incomplete  and  in  fact  on  continuing  the  substitu¬ 
tion,  contributions  from  the  displacement  equation  may  result  in  the 
cancellation  of  that  in  (41).  At  this  time  the  complete  analytical  asymptotic 
analysis  of  (39), (40)  is  a  subject  of  ongoing  research,  one  on  which  we  would 
welcome  the  assistance  of  other  investigators. 

In  actuality  the  unknown  pressure  p  is  not  the  quantity  of  greatest 
physical  import  but  rather  the  energy  release  rate  that  results  from  it. 

In  order  to  obtain  this  crack  driving  energy  G  from  p  we  proceed  as  follows. 

We  represent  G  at  a  given  depth  z  with  Sanders'  integral  Z3&1  in  an  xy-plane, 
take  advantage  of  the  symmetry  of  our  configuration  to  restrict  the  contour 
of  integration  to  one  lying  in  the  half-plane  y  >  0,  then  select  the  contour 
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2  2  2 

as  a  quarter-circle  (x  +  y  -  R  t  x  >  0,  y  >  0)  together  with  its  radial  ray 
on  the  y  axis  (x  =  0,  0  <  y  <  R) .  Now,  provided  we  include  the  plane  strain 
solution  for  a  crack  opened  under  constant  pressure  together  with  a  hydrostatic 
tensile  field  for  an  uncracked  half-space  (a  field  which  leaves  G  itself 
unaltered),  integrating  by  parts  and  invoking  equilibrium  results  in  the 
contribution  from  the  circular  portion  of  the  contour  vanishing  as  the  radius 
of  the  quarter-circle  tends  to  infinity  (R  ®),  leaving 


on  x  *  0.  To  obtain  the  terms  in  the  integral  of  (A3)  we  take  the  plane 
strain  and  hydrostatic  tensile  fields  as  well  as  the  stresses  and  displace- 

f 

ments  obtained  by  integrating  the  pressure  distributions  p,  p  multiplied  by 
the  appropriate  Boussinesq  kernels  over  the  surfaces  9 ^H,  9 respectively. 
We  can  then  use  (43)  to  evaluate  the  variations  in  G  along  the  crack  front 
when  we  assume  a  given  crack-front  shape;  for  a  crack  front  to  be  critical 
though  with  a  constant  G  we  must  solve  the  free-boundary  problem  described 
earlier.  The  adaptation  of  equations  (39), (40)  to  accommodate  such  a  crack- 
front  search  is  straightforward  and  details  may  be  found  in  C28l. 

The  integral  equations  for  assessing  the  effects  of  presumed  residual 
compressive  stresses  can  also  be  easily  derived,  essentially  by  following 
the  same  procedure  as  for  the  basic  infinite  3-D  problem.  That  is,  after 
subtracting  the  plane  strain  problem  to  ensure  convergence,  superposition 
of  point  loads  to  satisfy  the  outstanding  boundary  conditions  yields 
Integral  equations  identical  to  the  basic  problem  except  equation  (39)  is 
altered  by  adding  a of  (19)  to  its  right-hand  side. 

For  any  of  these  preceding  integral  equation  sets  our  approach  to  their 
numerical  solution  is  a  simple  one  -  a  collocation  procedure.  In  sum  this 
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entails  taking  due  account  of  the  expected  singular  behavior  of  the  unknown, 
transforming  it  so  as  to  admit  representation  as  a  piecewise  constant  on 
subregions  of  the  boundary,  and  thereafter  specializing  the  integral  equations 
to  hold  at  the  centroids  of  these  discretization  regions  to  generate  a 
square  linear  algebraic  system  of  equations  with  integral  coefficients  for 
the  determination  of  the  unknown. 

To  implement  the  procedure  we  first  nondimensionalize  the  integral 
equations  by  dividing  the  coordinates  x,z  and  variables  of  integration 
o,$,x,'P  by  the  semi-crack  length  a  (or  lim  a(z)  in  the  event  of  a  free 

Z-V«o 

boundary  problem),  and  the  stress  aP  (and  or  for  a  residual  stress  problem) 

z  y 

and  the  unknown  pressure  p  by  the  crack  opening  pressure  o°,  denoting 
all  the  quantities  so  formed  by  their  original  symbols  with  a  bar  atop. 

To  regularize  the  new  unknown  p  it  is  only  necessary  that  we  remove 
an  overestimate  of  its  singular  character.  Recalling  our  earlier  discussion 
in  the  introduction  it  would  seem  likely  that  this  could  be  achieved  in 
the  basic  problem  by  introducing  a  square-root  singularity  and  accordingly 
we  do  this  here  by  setting 

p  *  (x  >  1),  (44) 

Ax  -1) 


where  q  is  the  new  "regular"  unknown  and  the  particular  way  in  which  the 

square-root  singularity  is  introduced  in  (44)  is  in  anticipation  of  a 

subsequent  change  of  variable  to  remove  it.  In  addition,  from  our  previous 

2  2  2  2 

remark  that  the  decay  of  p  is  0(l/(x  +  z  ))  as  (x  +  z  )■*<*>  we  require, 

away  from  the  origin,  a  further  transformation  on  p  to  reflect  this  decay 
and  ensure  convergence  and  so  take 


JL£. 


(x2  +  S2)Ax2-l) 


-2  -2 
(x  +  z  »  1) 


■A/j-'A.Vl.V-.f- 


(43) 


Of  course  we  cannot  guarantee  that  q  of  (44),  (45)  is  indeed  continuous 
and  should  a  stronger  singularity  participate  at  the  crack  vertex  it  would 
not  be;  though  we  expect  this  is  not  the  case,  the  vindication  of  the 
regularization  in  (44) ,  (45)  ultimately  rests  with  the  numerical  results 
found  using  it.  In  contrast  we  do  expect  q  not  to  be  continuously  differ¬ 
entiable  at  the  crack  vertex  so  that  any  discrete  representation  other  than 
a  piecewise  constant  would  not  appear  warranted.  With  this  representation, 
the  discretization  of  the  surfaces  320.*  3^0.  on  which  q  is  represented 
proceeds  as  follows. 

Unlike  most  boundary  Integral  methods  the  surfaces  on  which  our 

equations  hold  are  not  of  finite  extent.  As  a  result,  the  initial  division 

of  the  boundary,  82^  U  3 32.*  *s  *nto  finite  regions  S  and  infinite  regions  I. 

For  our  preliminary  discretization  grid  or  map,  the  limits  of  the  finite 

regions  are  chosen  to  enclose  the  area  having  more  than  one  percent  variation 

from  the  normal  stress  component  o ^  of  the  plane  strain  solution  as 

estimated  by  the  previous  numerical  analyses  of  Cruse  and  Van  Buren  C243  and 

Tan  and  Fenner  C25].  Inside  the  border  of  this  region  the  grid  is  refined 

near  the  free  surface  and  crack  front  to  capture  the  expected  behavior  in 

the  unknown  q.  A  solution  to  this  system  is  found  for  an  initial  grid 

and  the  element-to-element  variation  in  q  used  as  an  indication  of  how 

the  grid  can  be  rearranged/extended  so  as  to  smooth  and  model  variations. 

This  partially  optimized  mesh  is  then  systematically  refined  to  study 

convergence  on  a  sequence  of  grids  comprised  of  a  coarse  (31  boundary 

elements),  a  medium  (109),  and  a  fine  (421).  To  indicate  the  maximum 

* 

resolution  the  fine  grid  is  shown  in  Fig.  6. 

The  only  remaining  numerical  issue  is  the  quadrature  to  be  used  to 

*Refer  Burton  C373  for  a  fuller  account  of  the  discretization  process. 


determine  the  matrix  coefficients  in  the  linear  set  of  equations.  This 
quadrature  must  recognize  the  singular  or  near-singular  behavior  of  the 
Boussinesq  response  kernels  as  well  as  the  singularity  introduced  through  (44) , 
(45),  the  latter  becoming  smooth  away  from  the  crack  front  (x  >  1) .  Thus,  to 
simplify  the  integration,  no  transformation  is  used  away  from  x  -*■  1+.  In  Fig.  6 
subdomains  of  S»  (i  =  1-3),  and  I,  (i  =  1-5),  are  shown  and  represent 
to  scale  the  sections  where  the  following  respective  transformations  are  used: 

1,  x//(x2  -1),  1;  1/z2,  x/(z2/(x2-l)) ,  1/z2,  l/(x2  +  z2) ,  1/x2.  A  consequence 
of  selectively  transforming  the  unknown  is  that  the  quadrature  on  S^, 
which  together  comprise  the  greater  part  of  S,  merely  entails  determining 
the  response  of  an  elastic  half-space  to  a  constant  pressure  applied  over 
a  rectangle  on  its  surface.  Hence  Love's  solution  C383  can  be  used  to 
perform  this  integration  exactly.  The  quadrature  on  the  remainder  of  S, 
the  surface  wherein  the  transformation  (44)  is  applied,  is  undertaken 
using  an  adaptation  of  the  first  mean  value  theorem  which  recognizes  the 
singularity  introduced  via  (44) .  The  quadrature  on  the  infinite  surfaces 
is  less  important  since  the  grid  itself  is  designed  so  that  the  unknown 
values  associated  with  these  regions  are  practically  negligible,  i.e.t 
so  that  the  key  results  should  be  left  virtually  unchanged  if  we  set  these 
integrals  to  zero.  With  this  in  mind,  the  only  singularity  directly 
handled  is  that  occurring  on  as  a  result  of  the  transformation  there, 
the  approach  being  the  advertised  change  of  variable,  specifically 
x  *»  /(x  -1).  Thereafter  all  of  the  infinite  surfaces  are  raappt  i  <to 
finite  ones  and  a  mid-point  rule  employed  -  a  rule  generally  consistent 
with  the  approximation  underlying  our  discretization.  For  those  few 
Instances  when  the  untreated  inverse-distance  singularity  in  the  Boussinesq 
displacement  kernel  makes  its  present  felt,  the  quadrature  is  confirmed 
simply  by  applying  the  mid-point  rule  on  successively  smaller  quadrature 


intervals. 


With  all  of  these  aspects  of  the  procedure  in  place  the  numerical 

solution  becomes  routine  since,  although  the  coefficient  matrix  is  neither 

symmetric  nor  sparse,  it  is  well-conditioned  and  quite  amenable  to  reduction 

by  standard  Gaussian  elimination.  While  the  computational  effort  required 

** 

is  not  great  for  the  resolution  gained,  the  results  found  directly  are 
limited  to  a  single  stress  component  in  the  crack  plane  and  even  for  this 
one  quantity  cannot  reasonably  be  expected  to  quantify  its  singular  character 
very  accurately.  Thus  it  is  in  part  in  order  to  overcome  these  short¬ 
comings,  as  well  as  to  provide  a  verifying  comparison,  we  next  turn  to  our 
finite  element  analysis. 

To  facilitate  grid  gradation  our  finite  element  approach  employs  two 
types  of  regular  or  host  elements  which  are  mutually  compatible:  an 
eight-node  isoparametric  hexahedron  or  bricky  and  a  six-node  isoparametric 
pentahedron  or  wedge  formed  by  coalescing  two  node-pairs  in  the  brick 
element  (Figs.  7A,  B) .  In  natural  coordinates  (£,n,C)  the  displacements 
U  are  represented  by  the  following  interpolation  functions 


u  =  c1  +  c2£+c3n  + V  +  +  c6n?  +  c +  cg£nc> 

u  =  ci  +  c2C+  c3?  +  c4(i+0n  +  c5&  +  c6(i+Onc> 


(46) 


Again  more  detail  can  be  found  in  Burton  E37D. 


**The  central  processing  times  on  the  time-sharing  DEC-20  system  at  Carnegie- 
Mellon  being  of  the  order  of  1  minute,  20  minutes  and  3  hours  for  the  coarse 
medium  and  fine  grids  respectively,  with  matrix  reduction  consuming  only  ~ 8 
minutes  of  the  fine  grid  time. 


where  (i  =  1-8)  are  vectors  of  constants,  readily  expressed  in  terms 

* 

of  nodal  values  of  U.  Being  isoparametric  the  respective  coordinate 
transformations  share  the  same  representations  as  in  (46). 

The  class  of  crack  problems  at  hand  also  requires  two  basic  types  of 
singular  element,  one  with  a  line  singularity  for  the  crack  front,  the 
other  with  a  point  singularity  for  the  crack  vertex.  We  construct  these 
singular  elements  after  Tracey  C39]  by  inserting  mid-side  nodes  in  a 
regular  element  and  letting  the  additional  degrees  of  freedom  model  the 
desired  singular  behavior.  The  advantages  of  this  method  of  development 
are  that  a  singular  element  so  formed  retains  its  constant  stress  and  strain 
capabilities  while  maintaining  full  interelement  displacement  continuity 
with  non-singular  elements.  On  the  negative  side,  the  coordinate  trans¬ 
formations  remain  those  of  the  original  regular  element  so  that  a  singular 
element's  edges  must  be  kept  straight  with  the  extra  nodes  fixed  at  their 
midpoints  in  real  (x,y,z)  space,  a  subparametric  element  therefore. 
Nevertheless,  provided  we  abide  by  these  restrictions,  such  singular  elements 
would  seem  to  be  among  the  best  for  present  purposes. 

In  detail  we  proceed  as  follows.  We  form  a  ten-node  pentahedron 
from  the  regular  pentahedron  (Fig.  7c)  and  use  the  extra  degrees  of  freedom 
to  represent  the  inverse  square-root  singularity  expected  at  the  crack 
front  away  from  the  upper  surface.  Hence  for  this  element 

u  =  cx  +  c2z+  c3c  +  c4(l+On  +  c5£  +  c6(l+0nc 


+  c7(l-O/(l+0  +  c8(l+oAl+0  +  c9(l-c)nAl+0  +  c10d+OnAl+£) . 

We  then  collapse  the  two  nodes  on  the  line  singularity  into  one  to  provide 


To  avoid  the  introduction  of  a  number  of  distinguishing  notations,  the 
one  nomenclature  is  used  in  (46)  and  future  interpolation  functions 
even  though  such  C  are  in  general  different. 


a  nine-node  tetrahedron  or  pyramid  (Fig.  7d)  and  free  up  the  singular 
nature  by  taking 


u  -  C1  +  c25+  c3(l+On  +  c4(l+0?  +  c5(l+On? 

+  c6(l+01_X  +  c7n(l+01_X  +  Cgtd+O1^  +  c9nc(i+01_X, 


(48) 


with  X  continuing  its  role  as  a  singularity  exponent.  The  corresponding 
singular  stress  behavior  generated  by  (47) ,  (48)  is  then 

CT  =  0(  if  Jr)  as  r  0,  0"  =  0(l/pX)  as  p  0,  (49) 

* 

r,  p  being  as  in  our  initial  coordinate  system  (Fig.  1) . 

To  determine  the  singularity  exponent  we  adopt  the  procedure  in  Swedlow 
C15J  wherein  X  becomes  an  additional  parameter  to  aid  in  minimizing  potential 
energy.  Of  course  this  minimization  seeks  the  lowest  potential  energy 
possible  using  the  limited  representations  available  in  the  singular  element 
and  thus  may  elect  to  model  some  part  of  the  regular  local  fields  that  is  not 
perfectly  captured  by  the  constant  strain  terms  with  a  little  of  the  X -sing¬ 
ularity  fields:  given  the  different  natures  of  the  two  fields*  such  eigen¬ 
function  smearing  is  not  seen  as  a  major  impediment  to  the  evaluation  of  X 
using  minimization  but  it  does  serve  as  a  warning  not  to  expect  tremendous 
resolution  in  such  determinations  in  complete  problems  that  have  quite  a 
number  of  different  local  fields  excited. 

The  quadrature  required  for  the  minimization  is  performed  using  Gauss 
rules,  with  singular  Gauss  rules  for  the  singular  elements.  The  actual 
schemes  used  are  extensions  to  3-D  of  the  two-dimensional  singular  rule 


To  ensure  compatibility  between  the  displacements  of  (47),  (48)  we  actually 
need  a  further  element  with  a  line  singularity,  a  transition  element  wherein 

the  /(1+0  in  the  Cg,  terms  is  replaced  by  (1+0 1  X  -refer  Solecki  [403 
for  greater  detail. 


used  by  Stern  C4ll.  These  extensions  are  simple  in  concept  but  somewhat 
algebraically  complex  because  of  the  higher-order  dimension  -  a  full 
description  is  given  in  Solecki  and  Swedlow  £423. 

With  the  representations,  minimization  procedure,  and  quadrature  in 
place  it  remains  to  apply  our  finite  element  method  to  the  class  of 
problems  of  concern  here.  As  a  result  of  some  numerical  experimentation, 
the  element  map  selected  for  the  baseline  problem  is  as  shown  in  Fig.  8  and 
has  375  elements  including  49  line-type  singular  elements  and  7  point-type, 

1587  degrees  of  freedom  all  told,  and  similar  resolution  to  the  medium 

* 

grid  in  the  integral  equation  analysis.  Perturbations  of  this  map 
also  serve  to  investigate  curved  critical  crack  profiles.  As  is  character¬ 
istic  of  FEM,  the  attendant  calculations  provide  a  wealth  of  information 
concerning  the  field  quantities  throughout  the  plate  from  which  it  is 
straightforward  to  calculate  energy  release  rates  using  a  path  independent 
integral  such  as  the  J  integral  (Rice  £433) .  Alternatively  these  key 
quantities  can  be  calculated  using  the  more  computationally  convenient 
device  of  virtual  crack  extension  (as  in  Parks  £443) ,  provided  care  is 
exercised  in  controlling  truncation  errors.  The  results  of  such  calculations 
together  with  those  from  the  integral  equation  analysis  are  considered  next. 

4.  Verification:  LEFM  implications  of  results 

Here  we  discuss  those  calculations  which  serve  to  validate  our  twin  analyses 
starting  with  checks  on  the  integral  equation  treatment,  then  comparing 
Integral  equation  results  with  finite  element,  and  finally  looking  at 

* 

CPU  times  for  this  grid  are  of  the  order  of  1  hour  compared  to  20  minutes 
for  the  medium  grid  in  the  integral  equation  approach  on  the  same  computing 

system. 


other  means  of  confirming  the  FEM  analysis.  With  the  resolution  of  the 
two  approaches  calibrated  we  close  the  section  by  considering  selected 
results  which  address  the  issues  raised  earlier  and  are  obtained  from  which¬ 
ever  method  is  deemed  most  advantageous  computationally. 

A  first  point  of  corroboration  of  our  integral  equation  approach  is 
found  in  the  magnitudes  of  the  three-dimensional  correction  away  from  the 
crack,  namely  the  values  of  the  unknown  p  on  the  infinite  surfaces  I^Ci^l-S). 
These  p-values  at  the  collocation  points  on  are  uniformly  less  than  2% 
of  o°,  thereby  vindicating  our  choice  of  the  extent  of  the  finite  regions 
in  the  discretization  and  supporting  the  justification  advanced  in  Section  3 
for  performing  the  quadrature  in  the  infinite  regions  with  a  simpler  but 
less  refined  scheme  than  in  the  finite.  A  second  point  of  corroboration 
lies  in  the  estimates  of  the  "regular"  unknown  q  which  are  indeed  numerically 
consistent  with  a  continuous  unknown  and,  in  particular,  exhibit  no  evidence 
of  a  stronger  singularity  at  the  crack  vertex  than  that  removed  by  the 
transformations  (44),  (45),  if  anything  indicating  the  opposite.  These 
numerical  results  therefore  add  credence  to  our  expectation  concerning 
singular  behavior  expressed  at  the  outset  and  also  suggest  that  the  numerical 
values  of  the  unknown  determined  via  the  integral  equations  should  converge 
well. 

To  examine  convergence  we  consider  the  unknown  pressure  p  by  itself 
rather  than  the  complete  stress  solution  including  the  plane  strain  response, 
since  it  is  for  this  quantity  alone  that  convergence  is  of  primary  import¬ 
ance.  Further  we  confine  attention  to  the  region  nearest  the  crack  front 
since  this  is  the  region  wherein  p's  role  is  most  significant.  Fig.  9 
shows  coarse,  medium  and  fine  grid  determinations  of  dimensionless  pressure 
p/a®  as  a  function  of  depth  z/a  at  the  closest  distance  from  the  crack  front 


COARSE 

MEDIUM 

FINE 


iation  unknown  p  of 
;k  (x/a  s  1.05,  v  = 


0.3). 
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common  to  all  three  grids,  x/a  =  1.05,  and  for  Poisson's  ratio  v  =  0.3. 

The  curves  in  Fig.  9  are  consistent  with  a  numerically  convergent  procedure, 
the  results  for  the  medium  and  fine  grids  being  almost  indistinguishable 
on  the  scale  that  space  permits  in  the  figure. 

As  the  finite  element  analysis  is  less  amenable,  in  terms  of  compu¬ 
tational  effort,  to  the  further  refinement  necessary  to  establish  conver¬ 
gence,  we  choose  instead  to  infer  this  quality  by  comparison  with 

integral  equation  results,  again  focusing  on  the  critical  region  near 
* 

the  crack  front.  We  do  this  by  taking  the  complete  normal  out-of-plane 
stress  in  the  crack  plane  from  the  integral  equation  results  in  conjunction  with 
the  same  stress  component  from  the  plane  strain  solution  for  a  crack  under  far- 
field  applied  tension ,  and  comparing  the  stress  so  assembled  with  the  correspond¬ 
ing  response  computed  using  finite  elements  on  plates  whose  plate  thicknesses  are 
Increased  until  the  mid-plate  answers  are  within  3%  of  the  integral  equation 
results  at  the  same  location  (the  actual  thickness  thus  found  being  h/a  =  2) . 

The  results  as  functions  of  depth  z/a  at  the  two  lines  of  centroids 
closest  to  the  crack  front  available  from  the  FEM,  x/a  =1.05,  1.14,  and 
for  v  =  0.3,  are  shown  in  Fig.  10  and  represent  good  agreement  given  the 
gradients  present  (c/.  for  example,  the  different  numerical  results  for 
the  same  configuration  in  Raju  and  Newman  El93>  P*  34).  Accordingly  it 
would  appear  that  the  finite  element  results  from  the  grid  of  Fig.  8 
are  convergent . 

Other  ratifications  of  the  FEM  analysis  are  the  perfect  response 
to  a  patch  test,  so  confirming  the  interelement  compatibility  claimed, 
and  the  evaluation  of  the  energy  release  rate  in  a  two-dimensional  test 

* 

In  this  connection  observe  that  halving  element  sides  to  form  a  more 
refined  grid,  the  procedure  in  essence  adopted  for  the  integral  equations, 
would  lead  to  an  FEM  map  of  around  3000  elements  and  13000  degrees  of 
freedom. 


problem  to  within  0.1%  using  the  element  arrangement  of  Fig.  8  and  the 
J  integral.  Though  not  strictly  a  verification  because  of  the  different 
physics  actually  present  and  that  theoretically  assumed,  an  experimental 
comparison  with  the  holographic  interferometry  study  of  Limtragool  C45U 
affords  an  additional  appraisal  of  the  resolution  of  our  FEM  procedure. 

Fig.  11  displays  w,  the  displacement  in  the  direction  of  the  crack  front  at  the 
crack  vertex  (see  Fig.  1  for  local  geometry,  Limtragool  C45],  p.  63  for 

full  specimen  geometry),  as  a  function  of  y,  the  surface  distance  transverse 

3 

to  the  crack,  for  high-strength  Aluminum  Alloy  6061-T6  (u  =  26  x  10  MPa,  v  *  0. 

both  as  found  experimentally  and  via  an  analysis  of  the  actual  specimen  geometry 

using  our  finite  element  capability.  The  agreement  shown  is  not  as  good  as  for 

the  same  displacment  component's  variation  with  the  other  surface  distance  x 

(Limtragool  C45j  p.  89)  but  is  nontheless  remarkable.  Also  presented  in 

Fig.  11  is  an  approximate  two-dimensional  estimate  as  determined  from 

simple  integration  of  the  out-of-plane  strain  in  a  plane  stress  analysis: 

the  discrepancy  between  this  last  and  the  physical  response  emphasizes 

the  inherent  three-dimensionality  of  the  displacement  plotted  in  Fig.  11. 

With  our  two  approaches  checked  we  now  consider  the  implications 

of  their  results  within  the  context  of  linear  elastic  fracture  mechanics, 

beginning  with  the  initial  step  in  LEFM,  the  identification  of  singular 
* 

nature.  For  this  activity, only  the  finite  element  approach  can  be 
reasonably  used  and  the  results  found  using  it  for  the  adjustable  singular¬ 
ity  exponent  X  of  (48),  (49)  as  a  function  of  Poisson's  ratio  v  are  sum- 
arized  in  Table  1,  together  with  those  of  other  researchers. 

* 

More  comprehensive  results  than  those  reported  here  may  be  found  in 
Burton  C37D  and  Solecki  C46U  for  the  integral  equation  and  finite 
element  investigations  respectively. 


Figure  11.  Comparison  of  FEM  and  experimental  results  for  the 
normal  displacement  on  the  free  surface  (Limtragool  [45]). 


TABLE  1 


Singularity  exponent  at  the  crack  vertex  (X  in  <J  =  0( 1/p*)) 


Source 

v  =  0.0 

v  *  0.15 

v  =  0.3 

v  =  0.4 

Analytical  by 

Benthem  C5] 

0.500 

0.484 

0.452 

0.413 

Analytical  by 

Kawai  et  al.  C6] 

0.50 

0.48 

0.43 

0.37 

Finite  difference 
by  Benthem  CllD 

0.500 

0.484 

0.452 

0.414 

FEM  by  Bazant 
&  Estenssoro  Cl2] 

0.500 

0.484 

0.452 

0.413 

FEM  by  present 
treatment 

0.499 

0.485 

0.445 

0.370 

The  results  in  Table  1  from  our  analysis  are  for  a  plate  with  h/a  =  1, 
there  no  longer  being  a  need  to  simulate  the  infinite  half-space  and  a 
less  distorted  grid  being  preferred.  An  indication  of  a  probable  upper 
limit  on  the  resolution  of  our  results  is  given  by  the  X -value  for  v  =  0 
which  is  0.001  in  error  from  the  known  exact  result  of  1/2.  In  view 
of  their  likely  resolution  for  other  v  they  should  not  be  interpreted 
as  supporting  the  cited  X -values  of  Kawai  et  al.  C6l  over  those  of 
Benthem's  analysis  C5],  but  rather  as  being  in  general  agreement  with  both, 
no  evidence  of  the  stronger  singularity  of  C63  being  found  in  our  specific 

configuration.  In  contrast,  the  close  agreement  of  the  numerical  work  of 
Benthem  CllD  and  Bazant  and  Estenssoro  Cl2l  with  Benthem’s  analytical  values  C5D 
does  support  the  last  over  the  corresponding  ones  in  Kawai  et  al.  C63  by  virtue 
of  the  fact  that  Cll],  Cl2]  treat  problems  especially  tuned  to  X -determination 
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instead  of  the  global  problems  analyzed  here.  The  import  of  all  the 
results  in  Table  1  is  that  the  energy  release  rate  should  be  zero  right 
in  the  upper  surface  for  normal  crack  intersection  when  v  ^  0  and,  given 
a  continuous  variation  along  the  crack  front,  that  G  should  decay  near  the 
surface  in  this  instance.  We  investigate  this  proposition  next. 

As  argued  previously  these  energy  release  rates  are  to  be  determined  in 
planes  parallel  to  the  free  surface.  The  method  favored  as  being  the  most 
reliable  for  this  determination  is  by  means  of  path  independent  integrals  in  such 
planes.  Since  the  FEM  analysis  yields  the  information  needed  in  such 
evaluations  more  readily  than  the  integral  equation  approach,  we  use  the 
finite  element  results  to  this  end  here.  Specifically  we  again  take  h/ a  =  1 
and  use  the  J  integral  (after  Rice  C433)  as  well  as  virtual  crack  extension 
(after  Parks  C443)  to  compute  G,  the  former  serving  to  check  the  more 
computationally  convenient  latter.  In  instances  in  which  both  methods  of 
G-determination  are  applied  the  same  results  are  obtained  for  all  practical 
purposes,  justifying  the  use  of  the  second  method  for  intervening  results. 

These  results  as  functions  of  depth  for  various  Poisson's  ratios  are  given 
in  Fig.  12,  with  G  therein  being  normalized  by  its  plane  strain  value  G^. 
Mathematically  the  indications  are  that  G  should  be  identically  zero  in  the 
free  surface  for  v  t  0  but  that  it  can  be  non-zero  arbitrarily  close  to  the 
free  surface:  numerically  we  find  a  decay  of  G  as  the  free  surface  is  approached 
(z/a  -*•  0)  with  the  largest  drop  between  the  maximum  value  and  that  of  the 
closest  Gauss  point  to  the  free  surface  occurring,  as  expected,  when  v  takes 
on  the  maximum  value  considered  (0.4),  and  being  30%  of  Gp. 

The  configuration  analyzed  to  provide  the  results  in  Fig.  12  does  not 

Regards  comparing  the  present  finite  element  analysis  with  that  of  Bazant 
and  Estenssoro  p.23,  it  is  worth  noting  that  0.2  3  has  the  highly  comparable 
results  of  X  *  0.485  for  v  ■  0  and  X  =  0.380  for  v  =  0.4  on  the  fine  grid 
therein  prior  to  exploiting  convergence  to  extrapolate,  an  option  not  open 

to  us  here. 


coincide  with  any  of  the  standard  specimens  used  in  the  determination  of 
the  critical  energy  release  rate  Gc>  i.e. ,  in  fracture  toughness  testing 
(see  C47j) .  Nonetheless  it  is  somewhat  similar  in  that  it  entails  opening 
or  Mode  I  excitation  and  in  that  h/a  *  1  is  the  nominal  ratio  prescribed 
for  all  of  the  standard  test  pieces  in  C47D.  Accordingly  we  can  infer 
what  inaccuracies  are  introduced  in  standard  toughness  tests  as  a  result  of 
not  accounting  for  three-dimensional  variations  with  the  two-dimensional 
plane  strain  value.  From  Fig.  12  we  see  that,  if  indeed  the  crack  front  is 
perfectly  straight  as  assumed  in  the  standards  E47],  then  the  Gc  calculated 
from  a  two-dimensional  analysis  of  a  test  could  be  as  much  as  18%  less  than 
the  real  material  value  for  v  =  0.4  and  typically  would  range  from  4-12%  smaller 
(corresponding  conservative  errors  in  the  critical  value  of  the  stress 

* 

intensity  factor  K£  being  9%  for  high  v  with  a  normal  range  of  2-6  %) . 

Given  the  overall  scatter  usually  encountered  in  fracture  toughness  testing, 
these  sorts  of  error  are  not  out  of  line. 

In  the  cyclic  life  estimation  component  of  LEFM  the  import  of  the 
results  of  Fig.  12  are  less  clear.  Ignoring  three-dimensional  effects  and 
assuming  a  single  A/G(AK)  for  a  through  crack  in  a  plate  in  a  cyclic  life 
calculation  is  certainly  not  completely  correct,  but  since  the  test  data 
used  in  such  calculations  is  usually  analyzed  with  the  same  simplifying 
assumption  the  errors  introduced  by  the  pair  of  simplifications  would  appear 
to  compensate  in  large  part.  Of  course,  3-D  LEFM  implies  that  for  v  ^  0 
the  straight  profiles  associated  with  the  varying  Gs  of  Fig.  12  would  not 
propagate  as  such  but  rather  take  up  curved  profiles  with  constant  energy 
release  rates:  we  next  investigate  what  shapes  these  crack  fronts  must 
adopt  in  order  to  compensate  for  the  drop-offs  in  G  displayed  in  Fig.  12. 

* 

These  statements  presume  that  the  specimen  would  in  fact  fracture  when  the 
three-simensional  maximum  G  attained  the  critical  material  value  G  for  un¬ 
stable  crack  propagation,  then  taking  the  G-value  calculated  from  two-dimen¬ 
sional  analysis  at  that  load  as  Gc  gives  rise  to  the  lower  values  indicated. 


In  view  of  the  free  surface  reductions  in  energy  release  rate  it  would 
seem  reasonable  to  advance  the  crack  front  away  from  the  upper  surface  so  as  to 
undercut  the  material  ahead  of  the  crack  near  the  surface  thereby  increasing 
the  energy  release  rate  there  and  approaching  a  constant  G  profile.  Indeed 
the  physical  evidence  suggests  this  is  the  case  (as  in,  say.  Bell  and  Feeney 
C33D)  and  Bazant  and  Estenssoro  Cl2U  theoretically  determine  a  crack/surface 
intersection  angle  which  deviates  from  perpendicular  in  a  manner  consistent 
with  subsurface  crack  advance  and  which  strengthens  the  vertex  singularity 
to  the  point  of  enabling  a  nontrivial  G  to  exist  right  in  the  free  surface. 
Accordingly  we  look  for  critical  crack  fronts  of  this  type  here,  beginning 
our  search  by  treating  moderate  v-values  because  we  know  the  answer  for  v  =  0 
is  simply  the  straight  crack  with  normal  intersection.  In  this  way  the 
critical  profiles  of  Fig.  13  are  generated  for  succesively  larger  values  of 
Poisson's  ratio,  each  profile  shown  therein  having  less  than  a  2%  numerical 
variation  in  G  with  depth  z/a.  These  results  illustrate  that  only  a  minor 
amount  of  crack  front  curvature  is  needed  to  render  a  profile  critical  in 
an  LEFM  sense,  the  greatest  deviation  from  straight  occurring  for  the 
maximum  v  of  0.4  and  representing  an  advance  of  less  than  2%  of  the  semi¬ 
crack  length  a:  that  is,  a  30%  variation  in  G  is  overcome  by  an  order  of 
magnitude  less  change  in  crack-front  profile. 

The  implications  of  these  results  regarding  critical  G  determination 
or  fracture  toughness  testing  are  as  follows.  Logically  LEFM  requires 
that,  if  toughness  testing  is  to  take  into  account  three-dimensionality, 
the  governing  test  specifications  should  admit  limited  two-sided  deviations 
in  actual  specimen  crack  fronts  about  curved  profiles  that  are  everywhere 
critical  like  those  of  Fig.  13.  As  it  happens,  present  standards  C^7j 
already  do  this  since  they  permit  a  considerable  curvature  and,  moreover,  the 
allowable  scatter  is  such  that  the  most  curved  of  the  present  profiles  for 


v  =  0.4  is  placed  well  within  this  spread.  Indeed,  in  reflecting  large 
amounts  of  testing  experience,  the  most  recent  standards  C47H  allow  far 
greater  curvature  than  the  most  curved  profile  in  Fig.  13,  possibly  because 
of  the  presence  of  some  shallow,  surface,  residual,  compressive  stresses.  With¬ 
out  quantification  of  such  fields,  an  accurate  assessment  as  to  the  appropriate¬ 
ness  of  the  standards  in  this  regard  is  impossible,  but  subsequent  results 
reported  here  do  indicate  that  the  increased  curvature  permitted  in  C47]  is 
qualitatively  reasonable. 

Implications  for  cyclic  crack  growth  calculations  in  LEFM  are  less 
easily  arrived  at.  While  it  is  clear  that  only  minor  modifications  in 
crack  front  shape  are  necessary  to  produce  the  sort  of  profile  which  under¬ 
goes  self -similar  growth  -  an  implicit  assumption  in  much  of  current  two- 
dimensional  cyclic  life  estimation  -  it  is  not  so  obvious  how  to  simply 
adjust  two-dimensional  calculations  to  account  for  the  overall  three-dimensional 
increase  (here  about  10%)  in  G,  hence  A/G.  As  mentioned  earlier,  provided 
these  increases  are  comparable  in  the  calibrating  test  and  the  application, 
their  effects  should  cancel  to  some  extent.  However  in  instances  wherein 
there  are  significant  differences  between  the  two  situations  it  is  possible  that 
significant  losses  of  accuracy  could  accrue.  One  possible  type  of  fundamental 
difference  between  a  data  generating  test  specimen  and  the  configuration  to 
which  it  is  applied  is  in  relative  thickness  and  we  consider  this  effect 
next. 

Fig.  14  presents  FEM  results  for  energy  release  rates  normalized  with 
respect  to  the  plane  strain  value  as  functions  of  depth  for  v  =  0.3  and 
the  two  finite  thickness  extremes  treated  here,  the  fattest  being  that 
used  to  simulate  the  basic  infinite  geometry  and  having  h/a  =  2,  the  thinnest 
having  h/a  =  0.3.  The  variation  in  G  increases  as  thickness  decreases  with 
the  change  from  maximum  to  minimum  going  from  11%  of  G  when  h/a  =  2  to  20% 
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of  G  when  h/a  =  0.3.  This  trend  may  explain  in  part  the  disagreement  be- 

P 

tween  the  integral  equation  and  finite  element  results  in  Fig.  10  since  it 
is  possible  that  the  FEM  results  there  showed  a  somewhat  greater  free-surface 

reduction  due  to  finite  thickness  effects,  viz ,  the  possibility  that  h/a  =  2 
does  not  completely  simulate  an  infinite  half-space.  On  the  other  hand 
the  value  of  G  in  Fig.  14  at  maximum  depth  for  h/a  =  2  would  seem  to  support 
the  contention  that  this  thickness  of  plate  is  close  to  "infinite". 

While  the  drop-off  in  G  changes  in  Fig.  14,  the  average  value  is  in¬ 
creased  by  less  than  3%  as  plate  width  reduces  -  well  within  the  10%  increase 
predicted  in  two-dimensions  for  v  =  0.3  as  one  goes  from  plane  strain  to 
plane  stress  -  and  the  excess  of  G  over  the  two-dimensional  value  (G^) 
remains  nearly  constant.  Consequently  little  error  would  seem  to  be  intro¬ 
duced  in  two-dimensional  cyclic  life  calculations  from  ignoring  three- 
dimensional  elastic  effects  when  the  test  specimen  differs  markedly  from 
the  application  geometry  in  thickness. 

A  further  3-D  factor  possibly  affecting  cyclic  life  calculations  are 
residual  stresses  and  Fig.  15  shows  the  effect  on  the  energy  release  rates 
of  the  following  presumed  residual  stress  fields;  superficial  with 
decaying  to  1%  of  a°  at  z/a  =  0.05  (y  ■  92.1  in  (19)),  shallow  with  decay 
depth  z/a  =  0.3  (y  =  15.4),  medium  with  decay  depth  z/a  =  0.6  (y  =  7.7), 
and  deep  with  decay  depth  z/a  =1.2  (y  =  3.8).  These  results  are  readily 
calculated  using  the  integral  equation  approach  since  we  can  now  check  a 
local  stress-fitting  determination  of  G  -  the  easiest  for  the  integral  equation 
results  -  with  a  more  reliable  finite  element  calculation  using  a  path 

independent  integral.  This  check  for  the  instance  of  no  residual  stress 
gives  the  attendant  energy  release  rate  Gq  to  within  1%.  Normalizing  the 

other  energy  release  rates  in  the  presence  of  residual  stresses  Gf  by  Gq» 
we  see  in  Fig.  15  that  sizable  reductions  in  energy  release  rate  occur  in 
the  presence  of  residual  compressive  stresses,  both  at  the  surface  and  at 


Figure  15.  Presumed  compressive  residual  stresses  and  resultant 
relative  energy  release  rates  for  normal  intersection  (ur0.3). 
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depths  considerably  greater  than  the  decay  depths  of  the  corresponding 
residual  stress  field  -  clearly  a  three-dimensional  effect  and  not  something 
that  one  would  readily  estimate  from  any  approximate  two-dimensional  analysis. 
While  these  calculations  are  for  hypothetical ,  presumed,  residual,  stress  fields 
they  do  serve  notice  of  the  potentially  major  impact  of  residual  stresses  on 
cyclic  life  calculations  even  within  LEFM,  and  the  need  to  quantify  such 
fields  therefore  if  better  cyclic  life  estimates  are  to  be  determined. 

5.  Concluding  remarks 

The  consensus  emerging  from  the  literature  of  a  weaker  singularity  at  a 
normal  crack/surface  intersection  in  analyses  of  local  problems  is  further 
supported  by  the  present  investigation  of  some  global  crack/ surf ace  inter¬ 
section  problems.  A  degree  of  confidence  in  the  present  results  can  be  gained 
from  the  good  agreement  of  two  independent  numerical  analyses  and  the  per¬ 
formance  of  these  two  approaches  with  respect  to  a  number  of  checks.  The 
two  methods  also  provide  information  on  the  key  quantity  in  linear  elastic 
fracture  mechanics,  the  energy  release  rate,  for  a  variety  of  configurations  from 
which  the  following  conclusions  may  be  drawn. 

The  three-dimensional  implications  of  crack/surface  intersections  re¬ 
garding  fracture  toughness  testing  may  reasonably  be  disregarded  as  it  would 
appear  they  can  be  in  most  present  day  cyclic  life  calculations.  One 
exception  to  this  last,  however,  may  well  be  in  near-surface  residual  stress 
effects  and  the  idealized  pilot  study  presented  here  indicates  a  need  to 
understand  such  fields  better.  In  sum,  but  for  residual  stress  effects,  it  would 
seem  unlikely  that  any  results  of  great  practical  consequence  would  be 
furnished  from  the  further  study  of  elastic  crack/surface  intersection  pro¬ 
blems.  This  is  not  to  say  there  are  not  some  important  three-dimensional 
elastic  crack  problems  to  be  considered.  Indeed,  the  apparent  paradox  of 
cyclic  crack  growth  at  a/G  levels  for  which  the  maximum  Gs  attained  are  less 
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than  Gc  strongly  suggests  that  this  problem  is  three-dimensional,  with 
portions  of  the  crack  front  departing  from  a  simple  straight  profile  and 
as  a  result  experiencing  energy  release  rates  in  excess  of  Gc:  some 
credibility  is  lent  to  this  possiblity  by  the  present  investigation  since  the 
rearrangements  to  obtain  critical  crack  fronts  may  be  reinterpreted  as  indicating 
that  small  changes  in  crack  profiles  can  lead  to  far  larger  variations  in  G, 
whence  a/G.  Used  with  sufficient  care  the  sort  of  methods  developed  here 
may  be  able  to  investigate  this  suggestion  and,  if  successful,  could  lead 
to  a  physically  reasoned  explanation  of  fatigue  crack  growth  -  certainly  an 
understanding  which  is  potentially  of  major  practical  significance. 
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